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Considerable evidence suggests that anticipating sudden shifts from one state to another in 
bistable dynamical systems is a challenging task, examples include ecosystems, financial markets, 
complex diseases, etc. In this paper, we investigate the effects of additive, multiplicative and cross 
correlated stochastic perturbations on determining regime shifts in a bistable gene regulatory sys¬ 
tem, which gives rise to two distinct states of low and high concentrations of protein. We obtain the 
stationary probability density and mean first passage time of the system. We show that increasing 
additive(multiplicative) noise intensity induces regime shift from a low(high) to a high(low) pro¬ 
tein concentration state. However, an increase in cross correlation intensity always induces regime 
shifts from high to low protein concentration state. For both bifurcation (often called tipping point) 
and noise induced (called stochastic switching) regime shifts, we further explore the robustness of 
recently developed critical slowing down based early warning signal (EWS) indicators (e.g., rising 
variance and lag-1 autocorrelation) on our simulated time series data. We identify that using EWS 
indicators, prediction of an impending bifurcation induced regime shift is relatively easier than that 
of a noise induced regime shift in the considered system. Moreover, the success of EWS indicators 
also strongly depends upon the nature of noise. Our results establish the key fact that finding more 
robust indicator to forewarn regime shifts for a broader class complex natural systems is still in its 
infancy and demands extensive research. 


I. INTRODUCTION 

Theoretical as well as experimental studies have indi¬ 
cated that many complex systems under the influence 
of stochastic perturbations can undergo sudden “regime 
shifts” in which they abruptly shift to a contrasting state. 
Such shifts may occur in systems with alternative stable 
states [I| . Well studied examples of regime shifts include 
sudden collapse of ecosystems 0, the onset of collapse 
in mutualistic communities @ , abrupt climatic shifts Q , 
the crash of markets in global finance systemic fail¬ 
ures such as the epileptic seizures @ and even the erup¬ 
tive events in spreading fires Q . Despite rich advances in 
the theory of complex systems, understanding the mech¬ 
anism which trigger the onset of regime shifts in nature 
remains a challenge. 

There are mainly two types of regime shifts that can 
occur in systems with alternative stable states. One is 
critical transition which is associated with the bifurca¬ 
tion points (so called tipping points) Q and another is 
noise induced transition (also known as stochastic switch¬ 
ing) i- Much effort has been devoted in recent years 
in developing early warning signals of impending regime 
shifts between alternative stable states (lOl-Hq. Such 
early warning signals can have tremendous impact on 
managing natural systems by forewarning the systems 
under the threat of state shift, so that appropriate man¬ 
agement strategies can be initiated to prevent a catastro¬ 
phe. Recent progress in this direction suggests that a set 
of generic statistical indicators (e.g., increase in variance, 
autocorrelation) may forewarn an impending transition 
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in a wide range of complex systems A recent review 
on “dynamical disease” argued that an early detection of 
regime shifts even in the field of medical sciences, such as 
in cardiac arrh yth mia’s could be of some help to prevent 
sudden death [T^. However, these signals are mostly 
developed for the phenomenon of critical slowing down 
that arises in the vicinity of tipping points [isl [T^ . For 
purely noise induced transitions, such as we also exam¬ 
ine in this paper, there is an active debate about whether 
early warning signals can really be useful [l^ [l^ . 

In genetic regulatory systems, it is known that ran¬ 
dom cell-to-cell variations within a genetically identical 
population can lead to regime shifts between alternative 
stable states of gene expression (i.e., sudden transition 
in protein concentration) [ 2 O, [2l[. When the underly¬ 
ing genetic system contains regulatory positive feedback 
loop s, individual cells can exist in different steady states 
[^ . Some cells may, for example, live in the “off” ex¬ 
pression state of a particular gene, whereas others are in 
the “on” expression state [l^. These stochastic fluctua¬ 
tions in gene expression, commonly referred to as noise, 
have been proposed to cause transitions between these 
states [13, [ 2 ^ • A well known example of bistable gene 
expression with coupled stochastic transition is the in¬ 
duction of the lac operon in E. coli which results in the 
synthesis of protein ff-galactosidase required for breaking 
up sugar molecules and releasing energy to the cell [^ . 
Experimental study on ff-galactosidase shows that sud¬ 
den transition from unregulated (low level) to regulated 
(high level) P-galactosidase state of lac operon occurs at 
a critical point of an inducer concentration . A poten¬ 
tial application of early warning signals in gene expres¬ 
sion is to identify increased risk of sudden transitions in 
protein concentration and prevent complex disease onset 
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Recently the concept of regime shifts with associated 
early warning signals has been used in systems medicine 
M- The expectation is to foresee a sudden catastrophic 
shift in health condition which may result in a extreme 
transition to a disease state. Now a days it is believed 
that a detailed understanding of regime shifts in dis¬ 
ease onset will provide broad applications in the held 
of medicine [l^. Already detected sudden transitions in 
medical science include epileptic seizures , depression 
[13, pulmonary disease [3l|, diabetes mellitus [HI, etc. 
For the case of type 1 diabetes mellitus, [i-cells in the 
pancreas do not produce protein hormone insulin, which 
sometimes is a consequence of “switched off” state of 
HLA-encoding genes in the cell (“switched on” state of 
genes in /?-cell corresponds to activation of insulin pro¬ 
duction). This is due to the fact that genes carry the 
instructions that cells use for protein production. Hence, 
bistability via “switching on” or “off” states in “gene ex¬ 
pression” is also an important topic to study for regime 
shifts in systems medicine [^ . 

In this paper, we study a stochastic version of bistable 
gene regulatory positive feedback loop model [13 to ex¬ 
plore the robustness of early warning indicators of regime 
shifts, for both cases, the critical transition and noise in¬ 
duced transition. We investigate the effect of additive 
and multiplicative noise intensities, and cross correlation 
intensity between two noises on the model by calculating 
the probability density and potential function.We hnd 
that increasing the intensity of additive noise induces 
regime shifts from low to high protein concentration state 
and vice versa for an increase in multiplicative noise in¬ 
tensity. Whereas an increase in the cross correlation in¬ 
tensity from a negative to a positive value between the 
two noises induce regime shifts from high to low protein 
concentration state. We also compute mean first pas¬ 
sage time (MFPT) for escape over the potential barrier. 
We discuss how oue can regulate the production levels of 
protein. To this end, we apply the early warning signals 
of regime shifts [l3 on the simulated time series data of 
the stochastic model to examine how successfully one can 
forewarn regime shifts. Our work presents a novel frame¬ 
work for using early warning signals to identify regime 
shifts, and also their key limitations, in a gene regulatory 
system. Fiually, in the discussion section, we conclude 
the paper with a discussion of the main results together 
with the applicability and importance of our study. 


scription when bind to a responsive element (TF-RE) in 
the DNA sequence (see Fig. [1] for a schematic represen¬ 
tation). The transcription factor activator TF-A is re¬ 
ferred to as regulatory protein, is used to control genetic 
regulation and acts as a pathway mediating a cellular 
response to a stimulus. The TF-A constructs a homod¬ 
imer which then binds to TF-RE or DNA regulatory site. 
Gene incorporates a TE-RE and when homodimers bind 
to the TE-RE, transcription of TF-A is increased. Ho¬ 
modimer binding to responsive element is independent 
phosphorylation of dimers. However, transcription can 
be activated only by phosphorylated TF-A dimers. Some 
dimers phosphorylation will depend on activity of kinases 
and phosphatases which is controlled by external signals. 
Hence, this genetic circuit incorporates signal-activated 
transcription as well as positive feedback on TF-A syn¬ 
thesis. Let, X denotes the protein TF-A, D{D*) denotes 
the unbound (bound) state of the DNA promoter, then 
the equilibrium reactions can be written as: 



fei/y" 

D + nX ^ D*, 
k2 

D* D* + X, 

where R is the basal expression rate, V is cellular volume, 
X is degraded with a rate constant kdeg, n represents co- 
operativity in binding, binding(unbinding) of TF-A ho¬ 
modimer to DNA regulatory site with a rate constant ki 
(k 2 ) and kp is protein production rate of TF-A. Letting 
X = [A], d = [D] and d* = [D*], the rate equation of 
the evolution of concentration of TF-A can be written as 



dx ^ x” 

- = R -f CL - 

dt Kd + x” 


kdeg^i 


( 1 ) 


where a = kp{d + d*) is the maximum production(i.e., 
transcription) rate, Kd = k 2 /ki is the dissociation con¬ 
stant of TF-A from TF-RE. A more detail derivation of 
Eq. ([T]) is given in [^ . 


1 1 0 



II. A MODEL OF GENETIC REGULATION 
A. Deterministic description 

It is well known that autoactivating positive feedback 
loop is one of the simplest circuit motifs able to exhibit 
bistable states in gene regulatory process [ssl - l^ . In the 
circuit, a single gene encodes a transcriptional factor ac¬ 
tivator (TF-A), and the TF-A activates its own tran- 


FIG. 1. A schematic picture of genetic expression for an auto¬ 
activating positive feedback loop. The expression of gene 
leads to protein (TF-A) that after oligomerization binds to its 
own promoter (TF-RE), acting as an self activator. Degrada¬ 
tion of the protein is denoted by the slashed circle (0). 

We can rewrite the dimensionless version of Eq. o as 

dx x" 

— = r + a- -; 


( 2 ) 
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where x = , t = kdeat, o, = . and f = 

" ’ fedes vTfd’ 

■ The analyses in this work have been carried 

kdeg </Kd 

out using the above model with dimensionless variable 
and parameters. For certain parameter values, Eq. © 
has three equilibria, of which two are attractors and third 
is a saddle point intermediate between the two attractors 
(see Fig. [2). The saddle point works as a basin boundary 
between the two stable equilibrium points. A thorough 
analysis of the deterministic model is given in [s^ . Here¬ 
after, in Eq. m for sake of simplicity we use x, r and a 
in place of i, t, f and d respectively, and Eq. @ becomes: 


dx 

dt 


( 3 ) 


If n > 1 and 0 < r < l/3-\/3, then Eg. ( 151 ) exhibits 
bistability for a range of the parameter a . To avoid 
lengthy calculations, we now replace n = 2 in Eq. ©, in 
the rest of this paper. To model ©, we investigate the 



FIG. 2. Bifurcation diagram for the deterministic model 
(black) and the stochastic model (green). The parameter val¬ 
ues are r = 0.1 (for the deterministic model) and r — 0.1 with 
CTi = 0.1, (72 = 0.1 and A = 0.1 (for the stochastic model). 
Stable steady states are marked with continuous lines and 
unstable steady states are marked with dashed lines, respec¬ 
tively. The stationary probability distribution for the additive 
noise model, Eq. (0), for different values of the production 
rate a is shown in cyan-red-white scale (color bar, logarithmic 
scale). 

effects of additive and multiplicative noise on the alter¬ 
native steady states of the gene regulatory circuit. We 
make this attempt because in the presence of noise, a 
bistable system may trigger regime shifts. 


B. Stochastic description 


of noise the system will eventually converge to one of its 
two stable fixed points. In which fixed point it will con¬ 
verge depends upon the initial condition. However, the 
presence of noise in the system will cause fluctuations 
in the steady states [4l|, which may lead to switching 
between two different stable states [T^ or there can be 
sudden transition from one stable state to the other sta¬ 
ble state [l^. In order to introduce the multiplicative 
and additive noise terms in Eq. we consider the one- 
variable Langevin equation in the general form as: 

^ = f{x)+g{x)^{t)+r]{t), (4) 

where f{t) and ri{t) represent Gaussian white noise. 
These noise terms have the following statistical prop¬ 
erties: (^(t)) = = 0, (^(t)^(t')) = 2cri5(t - t'), 

{r]{t )ig{t') ) = 2 a 2 S{t - t'), and = {v{t)kit')) = 

2X^aia2 S{t — t'), where ui and (T2 measure the level 
of noise strengths of f{t) and ri{t) respectively, A is the 
cross correlation between them, and t and t' denote two 
different moments. 

(^) Maximum production rate (a) 
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FIG. 3. Example time series datasets exhibiting regime shifts: 
(a) Critical transition (associated with bifurcation) and (b) 
Stochastic switching (purely noise induced) for the considered 
stochastic model. See text for more details. 


It is well known that noise is inherent in any natural 
system. In this work, we consider that dynamics of the 
protein concentration (x) are affected by multiplieative 
and additive noises, similar to Hasty et al. [d^. These 
noise terms can induce sudden shifts in the concentra¬ 
tion of protein (x). In a bistable system, in the absence 


In Eq. dU the additive noise r]{t) alters the background 
protein production [d^ . It is also known that in gene ex¬ 
pression, transcription is a complex sequence of reactions 
[4^ . thus it is expected that this part of the gene regu¬ 
latory sequence is also to be affected by fluctuations of 
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many intrinsic or extrinsic parameters. This implies the 
fact that the transcription rate (a) can be considered as 
a random variable |4fll ]. To vary the transcription rate 
stochastically, we consider a —?► a + ^{t). Hence, with 
the aforementioned modification, Eq. ([3]) (for n = 2) be¬ 
comes: 


dx 

dt 


ax 


- X + 


1 + x^ 1 + x^ 

fix) + g{x)^it) + r){t) , 




(5) 


where, fix) = r + — x and g{x) = Hence, the 

noise fit) is multiplicative, as compared to the additive 
noise ri{t). 

In the next section we study the stochastic model 
Eq. dS]) through a combination of analytical and simu¬ 
lation techniques. Our main aim is to investigate the ef¬ 
fects of multiplicative noise intensity cri, additive noise 
intensity U 2 and cross correlation strength A between 
two noises on the regime shifts between the high and 
low protein concentration states. In the analytical tech¬ 
nique these effects are studied by calculating probability 
densities, potential functions and MFPT. The simulation 
technique is complementary to the analytical technique, 
showing how the dynamical properties are captured in 
our analytical results can be seen within individual re¬ 
alizations based on example parameter sets, and adding 
information about how observed protein concentrations 
are arranged in time in these examples. In simulations, 
we also produce time series exhibiting regime shifts that 
can be analyzed using the same techniques as could be 
applied to real time series data (see Figs. [31(a)- (b)). The 
example times series can be divided into two broader 
classes: (a) critical transition time series (Fig. jSKa)) and 
(b) purely noise induced transition time series (Fig.|3l(b)). 


III. RESULTS 

A. Fokker-Planck equation and stationary 
probability density fnnction 

We begin this section by writing down the Fokker- 
Planck equation for the evolution of probability density 
of the dynamical variable x [4lj. Let P{x,t) denotes the 
probability density, which is the probability that the pro¬ 
tein concentration attains the value x at time t. Then, 
the Fokker-Planck equation (FPE) of Pix, t) correspond¬ 
ing to Eq. m is given by (^ : 

f)P(X P 3 3 ‘^ 

-m = + ^[Bix)Pix,t)] , ( 6 ) 

where 

A{x) = fix) aig{x)g'ix) g'ix), and 

Bix) = CTi [gix)]'^ -I- 0-2 -I- 2Ai/(Ticr2 gix). 

The limit of Pix, t) as t —>■ oo yields the stationary prob¬ 
ability density function (SPDF) of x, which we denote as 
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FIG. 4. Potential landscapes demonstrating how changes in a 
system parameter can cause decrease in resilience of equilib¬ 
rium point. Recovery rates upon stochastic fluctuations are 
lower if the basin of attraction is small (d) than that of a 
larger basin of attraction (a). The effect of reduced resilience 
can be determined by stochastic fluctuations induced in a sys¬ 
tem state ((b) and (e)) as increased standard deviation (S.D.) 
and lag-1 autocorrelation ((c) and (f)). Data sets to plot this 
figure are generated from Eq. with r = 0.1, (Ji = 0.0, 
a 2 = 0.05: (a-c) a = 2.5 and (d-f) a — 1.9. 


Psix). The SPDF Psix), which is the stationary solution 
of the FPE in Eq. ([6|) is given by [i^ (see Appendix SI 
for more details) : 


Psix) = 


Ne 

Bix] 


■ exp 


^ Ajx') 
Bix') 

Nr. 


dx 


'Xi b(a;)] + 0 - 2 - 1 - 2 A./crio -2 gix) 


exp 


r fix') + aigix')g'ix') + Av^o-io -2 g'jx') 
J <Ti [gix')f + 0-2 + 2X,/afdd gix') 


( 7 ) 


where W is normalization constant. Equation © can 
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also be put in the form: 

, ( 8 ) 


where 


1 r 2 

4'{x) = 2 + ^2 + 2Ai/cricr2 g{x) 


f 


f{x')dx' 


O'! [s(2;0] + 0-2 + 2 A^(TifT 2 g{x') 


(9) 


is called stochastic potential of the system [4^. The 
potential function maps the equilibria of dynamical sys¬ 
tems and their basins of attraction, by analogy to a “en¬ 
ergy landscape” in which the system state tends to move 
“downhill”. Extending this concept to stochastic dynam¬ 
ical systems gives a probabilistic potential (j){x) that com¬ 
plements the SPDF Ps{x) in characterizing the asymp¬ 
totic behavior of the considered system. Next we calcu¬ 
late (/)(x) and Ps{x) for three different cases concerning 
the effects of multiplicative and additive noise on antici¬ 
pating regime shifts in gene expression. We also calculate 
bifurcation diagram (see Fig. [2|) with changing the maxi¬ 
mum transcription rate a for both the deterministic and 
stochastic model. Note that, there is an enlargement of 
the bistability region for the case of correlated noise. In 
Fig. [21 we depict the stationary probability distribution 
of the additive noise model for different values of a, which 
is shown in white-red-cyan colorbar. This gives an idea 
about how extrema of stationary probability distribution 
is changing with the parameter a. Now, using analyti¬ 
cal techniques we first determine the parameter space 
where the system persists bistability still in the presence 
of stochasticity. Then, for a specific set of parameter 
values, we simulate time series of the system and finally 
using EWS indicators Q, we will try to detect forthcom¬ 
ing regime shifts. 


1. Additive noise 

First, we present the effect of an additive external noise 
source on the regime shifts between high and low concen¬ 
trations level of protein in gene expression. Hence, only 
the additive noise term is present in Eq. m and we as¬ 
sume that CTi = 0. Bifurcation diagram and stationary 
solutions of the additive noise model are same as for the 
deterministic model [13 • For ui = 0, the FPE (l6|) can 
be written as: 

dPix tl 3 3‘^ 

= -^[/(x)P(x,t)] +a2^[P(x,t)], (10) 

and similar to Eq. m the stationary solution Ps{x) is 
written as: 

Ps{x) = , (11) 


where 

(j){x) = i In (72 - — [ f{x')dx' , (12) 

2 02 J 







FIG. 5. (a)-(c) The effect of maximum production rate a 
on the stochastic potential (j>{x), when only additive noise is 
present in the system: (a) a = 1.3, (b) a = 1.9, and (c) a = 2.7 
with r = 0.1 and (72 = 0.05. (d)-(f) The effect of additive 
noise strength CT 2 together with maximum production rate a 
on the SPDF Ps{x): (d) a = 1.3, (e) a = 1.9, and (f) a = 2.7 
with r = 0.1 and for three different noise strengths (72 = 0.05 
(black curve) , (72 = 0.1 (red curve) and (72 = 0.2 (blue curve). 


is the potential and Na is the normalization constant. 

Figures |D(a) and |D(d) show the effective potential at 
two different values of a with a non-zero additive noise 
intensity. The depth and width of a potential well can 
be related to its resilience. Resilience is defined as the 
ability of a system to recover to its original state upon 
a perturbation. For critical transitions as the system 
approaches to a bifurcation point (see the positions of 
the parameter a in Fig. [2] marked by circles) the recov¬ 
ery rate from perturbation decreases smoothly, known as 
critical slowing down and as a result the system loses its 
resilience. One important prediction is that the loss of 
resilience should lead to an increase in the standard de¬ 
viation and autocorrelation (see Figs. IDJc) andSKf)) [1|- 
Next, we have demonstrated how the probability distri¬ 
bution of a system can be related with its potential and 
acts as a mapping between potential and their resilience. 
We simultaneously try to see the effect of changes in the 
parameter a and noise intensity on regime shifts. 

The effect of changing the maximum production rate 
a on the stochastic potential (jj^x) and simultaneously 
changing a together with different additive noise inten¬ 
sity (72 on SPDF Ps{x) is depicted in Fig. O In the 
left panel of Fig. HI we show that by varying a from low 
(see Fig. HKa)) to high values (see Fig. HKc)), the system 
can pass from a monostable low protein concentration 
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state (Fig.[5ja)) through a region of bistability (Fig.[5|b)) 
(i.e. coexistence of low and high concentration states) to 
a monostable high concentration state (Fig. EJc)). For 
a = 1.9, the region of bistability is distinguished by the 
presence of two local minima in the potential <(>(a;) (see 
Fig-HKb))- As already stated, the depth and width of the 
potential well determine the resilience of the equilibrium 
point. It is evident from Fig. [Hb) that the high concen¬ 
tration protein state (x) has high resilience. Hence, in 
a stochastic environment it has a low risk to experience 
a regime shift and can sustain under large disturbance 
in comparison with the low concentration protein state. 
On the other hand in stochastic models, peaks of the 
SPDF Ps{x) correspond to attractors and troughs cor¬ 
respond to repellors. Moreover, an equilibrium point is 
more stable (high resilience) if the SPDF peak is large in 
comparison with another equilibrium point which is less 
stable (low resilience) as the SPDF peak is small. The 
effect of three different values of noise intensity (T 2 on 
SPDF Ps{x) is shown in Figs. [ni(d)-[nKf) for three differ¬ 
ent values of a. In case of bistability, for a fixed a = 1.9, 
if we slowly increase the noise strength a 2 y the system 
still has coexistence of low and high concentration pro¬ 
tein state for low value of (72 • As we further increase 
CT 2 , the system will slowly shift to a single state consists 
of high protein concentration only (see Figs. EJe)-!!])!)). 
This regime shift is expected as because for a = 1.9 the 
high protein concentration state has high resilience. 


2. Multiplicative noise 
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0.6 
0.5 
3 0.4 
'i.'o.s 
0.2 
0.1 


(d) 


.a,=0.15 
..c,=0.25 
..n,=0.65 








FIG. 6. (a)-(c) The effect of maximum production rate a on 
the stochastic potential <?!>(a;), when only multiplicative noise 
is present in the system: (a) a = 1.2, (b) a = 2.3, and (c) 
a = 3 with r = 0.1, cri = 0.65 and (T 2 = 0. (d)-(f) The effect 
of noise strength cri together with maximum production rate 
a on the SPDF Ps{x)\ (d) a — 1.2, (e) a = 2.3, and (f) a = 3 
with r = 0.1 and for three different noise strengths cri = 0.15 
(black curve), cri = 0.25 (red curve) and cri = 0.65 (blue 
curve). 


Next, we consider only the presence of multiplicative 
noise in the system, i.e., in Eq. ([S]) the strength of addi¬ 
tive noise a 2 = 0. For a 2 = 0, the FPE ([61) becomes: 


(13) 

where 

A{x) = f{x) -b (Tig{x)g'{x), and 
B{x) = ai [g{x)f . 

The stationary solution Ps{x) of Eq. (ESI) is: 

Ps{x) = (14) 


where 


(^(x) = 2 [^1 [9ix)Y 


f{x')dx' 
<xi [9{x')f 


(15) 


is the stochastic potential and Nm is the normalization 
constant. 

Similar to the case of additive noise in the previous sec¬ 
tion ISec. lIII AD) , here Fig.E]shows the effect of changing 
the parameter a and the noise intensity cri on the poten¬ 
tial function (t){x) and the SPDF Ps{x). In one hand. 


when we increase the parameter a with a fixed noise in¬ 
tensity (see Figs. [6Ka)-l6Kc)), the system passes through 
a low concentration state to a high concentration state 
via a bistable region (Fig. EKb)). On the other hand, 
with variations in a we plotted the SPDF for different 
values of the noise intensity cri (see Figs. [6Kd)-l6Kf)). For 
monostable regions the effect of noise intensity on Ps{x) 
is similar, with an increase in ci the SPDF peak is de¬ 
creasing. However, in the bistable region (Fig.[6Ke)), with 
an increase in cti, the SPDF peak at the low protein con¬ 
centration X is increasing and that of the high protein 
concentration x is decreasing. This actually supports the 
fact that when there is bistability in the system, the low 
concentration state is more resilient (more stable under 
perturbative condition) than that of the high concentra¬ 
tion state, which is also evident from the Fig. E^d). 


3. Additive and multiplicative noise with correlation 

Herein, we consider that both the additive and mul¬ 
tiplicative noise git) and ^(t) are present in the system 
Eq. (El). Moreover, both these noise are statistically cor¬ 
related with a cross correlation parameter A. The correla¬ 
tion can arise from the regulation of feedback mechanism, 
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FIG. 7. (a) The potential (?i(a:) and (b) the SPDF Pb{x) for 
three different values of A : —0.25, 0 and +0.25. The other 
parameter values are r = 0.1, a = 1.9, cri = 0.4 and (T 2 = 0.9. 


i.e., in the presence of noise the protein concentration x 
is chemically coupled to the transcription rate a. Here 
we discuss three cases corresponding to the degree of cor¬ 
relation A, which is considered negative, zero or positive. 
In this section, we are mainly interested in to investigate 
how the correlation between ^{t) and r](t) triggers regime 
shifts from high to low protein concentration state and 
vice versa. We also showed that by changing cross cor¬ 
relation parameter A one can regulate the production of 
protein levels . 

The expressions of the SPDF Ps{x) and the poten¬ 
tial (l){x) are given in Eq. ([T]) and Eq. ([9]), respectively. 
Figure [7] shows the results for (j){x) and Ps{x) for three 
different strengths of correlation between additive and 
multiplicative noise. The system parameters r and a are 
chosen in such a way that the deterministic system re¬ 
tains bistability. As we vary the correlation parameter 
A from a negative value to a positive value, it is evident 
from Fig. [7l[a) that the resilience of right potential well 
is more sensitive with the change in A than that of left 
potential well. In fact, negative correlation (A < 0) will 
increase the stability of the high concentration of protein 
and positive correlation (A > 0) will increase the stabil¬ 
ity of the low concentration of protein (see the SPDF in 
Fig. EKb)). Hence, an increase in the cross correlation 
intensity always induces regime shifts from high to low 
protein concentration state. 


B. Mean first-passage time 

In biological systems, where noise plays an important 
role in determining the dynamics, it is of general interest 
to calculate the robustness of the system steady state un¬ 
der perturbative condition. The robustness of a steady 
state can be quantified by the mean first-passage time 
(MFPT) in noise driven systems [13,113 • In the bistable 
potential of Eq. (O, let xf/x^ {xf < x^) be the two 
steady states corresponding to the low/high protein con¬ 
centrations separated by the potential barrier x™ (i.e., 
the unstable equilibrium point). An equilibrium point 
can exit from its potential well in the presence of noise. 
The exit time is a number which depends on the specific 
realization of the random process and is known as first 
passage time. When the first passage time is averaged 


over many realizations, we get the mean first passage 
time. 




FIG. 8 . (a) The logarithmic value of MFPT with an increase 
in a. The other parameter values are r = 0.1, cti = 0.2, 
CT 2 = 0.9 and A = 0.01. (b) The logarithmic value of MFPT 
with an increase in A. The parameter values are r = 0.1, 
a = 1.9, CTi = 0.4 and (T 2 = 0.9. 


In the context of anticipating regime shifts, MFPT pro¬ 
vides a very useful characterization of the time scale on 
which a critical transition is likely to happen. When 
the structure of a dynamical system is known, the early 
warning signal indicators of regime shifts should be mea¬ 
sured from the time series data of length shorter than the 
MFPT of the state in that particular regime. 



FIG. 9. Plot of the ratio 7 = of the MFPT with respect 

to the control parameter a. The other parameter values are 
r = 0.1, CTI = 0.4 and CT 2 = 0.9 and A = 0.1. 


It is clear that the basin of attraction of the state 
extends from x™ to + 00 , as it is in the right of xf. Let 
T(x) be the MFPT to state x™ starting at x > x™. 
Then T(x) satisfies the following ordinary differential 
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FIG. 10. Early warning signals for simulated time series data of the stochastic model in the case of critical slowing down. 
We calculated variance and autocorrelation within rolling window of half the length of the time series segment (a segment is 
indicated by shaded region): (a) Additive noise: parameter values are r = 0.1, (ti = 0 and (J 2 = 0.9; (b) Multiplicative noise: 
parameter values are r — 0.1, cri = 0.9 and cr 2 = 0; and (c) Correlated noise: parameter values are r = 0.1, cri = 0.4, <72 = 0.1 
and A = 0.5. See text for more details. 


equation [44 


A{x\ 


dT{x) 

dx 


2 dx^ 


= -l, 


(16) 


with boundary conditions T(a;^") = 0 and = 0 . 

Similarly one can calculate the MFPT to state a;™ for 
the basin of attraction of the state xf^ which extends 
from 0 to xY\ After solving Eq. m, one can obtain the 
MFPTs as M\: 


T{xf) = 2 


r«‘) = 2 



B{z) 

B{z) 


dz, and 


dz, 


where 


The effect of changing a and also the correlation pa¬ 
rameter A on the MFPT are shown in Fig. [5] (for the 
details of parameter values see the caption of Fig. El). 
We can observe that the MFPT T{xf) decreases and 
T[x^) increases with an increase in a (see Fig.EJa)) and 
the opposite trend is observed for an increase in A (see 
Fig. EKb)). These also imply the fact that an increase in 
a promotes the regime shift from left potential well to the 
right potential well and vice versa for an increase in A. 
As MFPT T{xf) approaches to zero at the bifurcation 
point, the system loses its stability at low protein con¬ 
centration. The asymmetry in the exit time is calculated 


by the ratio 7 = and it diverges at the bifurca- 

tion point where T{xf) approaches to zero (see Fig. E]) 


44| . The quantity 7 thus gives an early warning signal 


of regime shifts from low to high protein concentration 
state for an increase in a. 


'ip{x) = exp 



2A{x') 

B{x') 



C. Early warning signals 


with xo = 0 for the x®* —>■ xf transition and xq = x^^ Here we assess the robustness of early warning signals 

for the xf^ —> x^ transition. to forewarn upcoming shifts to an alternative regime by 
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FIG. 11. Early warning signals for simulated time series data of the stochastic model in the case of stochastic switching. Here 
also we calculated variance and autocorrelation within rolling window of half the length of the time series segment (a segment is 
indicated by shaded region): (a) Additive noise: parameter values are r = 0.1, a = 2.1, cri = 0 and CT 2 = 0.9; (b) Multiplicative 
noise: parameter values are r = 0.1, a = 2.3, a\ = 0.9 and (T 2 = 0 and (c) Correlated noise: parameter values are r = 0.1, 
a = 2.1, (Ti = 0.4, (T 2 = 0.65 and A = 0.5. See text for details. 


analyzing simulated time series data [Hill from the 
considered stochastic model of genetic regulation. The 
simulation approach is showing how the characteristics 
captured in our analytical calculations can also be seen 
within individual time series realizations based on specific 
parameter sets. Early warning signals are observable sta¬ 
tistical signatures that antecede some state shifts. These 
signals have mainly been derived from the phenomenon of 
critical slowing down (CSD) which arises in the vicinity 
of bifurcation wherein its dominant eigenvalue will cross 
zero Q. As the eigenvalue reduces to zero, the response 
of the system becomes very slow in recovering from per¬ 
turbations and certain statistical features, such as in¬ 
creased variance and lag-1 autocorrelation, are predicted 
to appear in the time series analysis as a result (reviewed 
in Scheffer et al. [13 )■ Although many regime shifts 
appear as sudden transitions to another state in natural 
systems, the actual fact that not all regime shifts are as¬ 
sociated with a bifurcation [3 13 . Hence not all regime 
shifts associated with different mechanisms are expected 
to exhibit CSD. In Dakos et al. , it is pointed out that 
CSD based early warning signals “are not a panacea for 
anticipating all types of regime shifts”. 

Below we present two mechanisms those are respon¬ 
sible for observable regime shifts in our simulated time 
series. One is associated with the saddle-node bifurca¬ 
tion and exhibits CSD and another is associated with 
stochastic switching (SS) (i.e., purely noise-induced tran¬ 


sition and not associated with bifurcation) and does not 
exhibit CSD. There has been a recent debate about the 
success of early warning signals for predicting stochastic 
switching induced regime shifts [31 • Keeping that in 
mind, we felt that it is worthwhile to present what oc¬ 
curs in our system. For our analyses, we consider three 
different stochastic time series, each for CSD and SS: 
first we consider the presence of additive noise only, sec¬ 
ond we consider the presence of multiplicative noise only 
and finally presence of both noises with correlation in 
the system. To obtain the time series, stochastic simu¬ 
lations were performed in MATLAB (R2011a) using the 
Euler-Maruyama method (47l | and a standard integration 
step-size of 0.001. 

In our simulated time series, we visually identify shifts 
between low to high protein concentration and vice versa 
for both the cases, CSD and SS (see Figs. fTOl and fTTIl. 
Then we took different time series segments (the gray 
shaded regions in Figs. [T0| and [TT]) of different lengths 
(keeping in mind that their time lengths should be less 
than their MFPTs) preceding a regime shift and ana¬ 
lyzed those time series for the presence of early warn¬ 
ing signals. The “Early Warning Signals Toolbox” 
(http://www.early-warning-signals.org/) is used to per¬ 
form the statistical analyses. To ensure stationarity in 
residuals, we used Gaussian detrending with bandwidth 
40, on the time series data before performing any sta¬ 
tistical analysis. Then using a moving window size of 
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half the length of the considered time series (i.e., 50% 
of the considered time series segment), we calculate the 
variance and autocorrelation in our state variable, x, as 
these two indicators are both very commonly applied to 
anticipate regime shifts. The autocorrelation at lag-1 is 
measured by the autocorrelation function (ACF): 

E[(x(t) - / 2 )(x(t + 1) - 
Pi — -5-: 

where x{t) is the value of the state variable at time t, and 
/i and are the mean and variance of xit). Variance is 
the second moment around the mean ^ and measured as: 

- Pf ^ 

' i=l 

where N is the number of observations within the consid¬ 
ered time frame. A concurrent increase in both of these 
indicators forewarn an impeding regime shift Q. 

In Table I, we summarize results of early warning sig¬ 
nals (see Figs. [Till and fTTl) for the considered two different 
cases of regime shifts with three subcases each. The sig¬ 
nals [variance (cr^) and autocorrelation (pi)] are indicated 
as “-I-” if there is a concurrent rise; indicated as ” 

if there is a rise, but some false pick in the signals before 
the regime shift due to the presence of large fluctuations 
in the data; and indicated as ” if looking at the sig¬ 
nals, it is not possible to forewarn an impending regime 
shift. The result in Table I shows that in the case of CSD 


TABLE I. Two mechanisms, CSD and SS that can produce 
regime shifts with and without bifurcation in the system and 
the outcome of early warning signals for regime shifts. 


Noise Type 

CSD SS 

Var(a^) ACF(pi) Var(a^) ACF(pi) 

Additive 

+ -b - - 

Multiplicative 

+ +/- + + 

Correlated: 
Additive and 
multiplicative 

+ +/- - 


with all three different types of stochasticity applied to 
model (O, the variance is robust and always successfully 
predict upcoming regime shifts for all the three types of 
noise. However, the lag-1 autocorrelation is only suc¬ 
cessful to predict regime shifts for the case of additive 
noise only and can’t predict very positively an upcoming 
regime shift for multiplicative and correlated noise. This 
establishes the fact that even for the case of CSD the 
variance is more robust than autocorrelation in detect¬ 
ing regime shifts in genetic regulation and it is line with 
the findings in [i^. For SS, both of the early warning 
signals are not very successful in detecting regime shifts. 
Using the EWS indicators it is not possible to detect up¬ 
coming regime shifts for additive and correlated noises 
(see Fig. [TT]) . However, in the presence of multiplicative 


noise, the early warning signals give positive results. This 
is indeed a nice result given the fact that early warning 
signals are not developed for noise induced transitions, 
but rather for CSD which is associated with bifurcations. 
Nevertheless, our results show that even for the case of 
CSD though the EWS somewhat gives positive results, 
still some regime shifts may not be triggered by a con¬ 
current rise in variance and autocorrelation. This is due 
to the fact these CSD indicators are reliable for specific 
systems and also have some key limitations [3. Like, 
there is false positive and false negative errors in data 
[^ . Some other limitations include size of the window 
which varies across the literature and how much data are 
required for analyzing EWS is still not clear. Moreover, 
variance is a robust indicator as compared to autocorre¬ 
lation in the case of CSD and this is due to the fact that 
autocorrelation is affected by the length of time series 
(see Appendix S2 for more details) and more sensitive to 
false alarms. In the case of stochastic switching, antici¬ 
pating regime shifts is more difficult due to sudden tran¬ 
sitions. The significance of these findings are discussed 
in the discussion section below. 


IV. DISCUSSION 

In this paper, we investigated a stochastic genetic reg¬ 
ulatory circuit using analytical techniques and numeri¬ 
cal simulation to anticipate regime shifts in protein con¬ 
centration. In this respect, we derived an approximate 
Fokker-Planck equation from the Langevin equation. We 
studied the effects of intensity of both the additive and 
multiplicative external noise and their correlation on the 
gene regulatory system. The present study suggests 
that the presence of additive noise in the system induces 
regime shifts from a low protein concentration to a high 
protein concentration state, whereas multiplicative noise 
induces regime shifts from a high to a low protein concen¬ 
tration state. Moreover, we also show that how a correla¬ 
tion between the additive and multiplicative noise is im¬ 
portant in determining regime shifts and hence can reg¬ 
ulate the production of protein levels. Furthermore, an 
increase in the cross correlation intensity from a negative 
to a positive value between the two noises induce regime 
shifts from high to low protein concentration state. Here 
we have also coinputed the robustness of steady states us¬ 
ing the MFPT [^. Our MFPT result uncovers the fact 
that an increase in MFPT of right potential well with the 
maximum production rate promotes regime shifts from 
left potential well to the right potential well and vice 
versa for an increase in the cross correlation intensity. 

In addition, we used a recently developed tool 
(http://www.early-warning-signals.org/) of early warn¬ 
ing signals to anticipate regime shifts in the considered 
gene regulatory system using simulated time series data 
for both the CSD and SS. Extensive literature on regime 
shifts is available for ecosystems, financial markets, cli¬ 
matic shifts, etc. However, to the best of our knowledge. 
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there is very less work available on anticipating regime 
shifts in many areas of developmental biology, specifically 
in genetic regulation [iilHi. Here for the first time, for a 
genetic regulatory system, we used the time series anal¬ 
ysis based EWS approach to predict upcoming regime 
shifts. As anticipated by some previous authors [9|, ll9|, 
we also find that the EWS of raising variance and au¬ 
tocorrelation are in general sensitive to false alarms and 
not always successful to reliably predict impending state 
shifts in our model. We found that the EWS are moder¬ 
ately more successful when we analyzed time series data 
in advance to a state shift in the case of CSD, whereas 
in the case of SS, it is specific to particular noise only. 
We observed some key limitations and statistical issues 
of EWS, such as false pick in data, size of the window 
and length of time series data. For our model, we also 
verified that rising variance appears as a robust indica¬ 
tor of CSD as compared to lag-1 autocorrelation. This 
is mainly due to the fact that autocorrelation is sensitive 
to the length of time series. For accurate estimation of 
autocorrelation, we need long and equidistant time se¬ 
ries data which are not always available for real systems. 
However, variance is insensitive to this effect and as a 
result, measuring autocorrelation as EWS may increase 
the possibility of false alarms. While testing EWS to our 
simulated time series data, we also observed that selec¬ 
tion of window size is another major problem for getting 
a positive signal of regime shifts. In the case of SS, an¬ 
ticipation of regime shifts is extremely difficult because a 
bifurcation point is neither approached nor crossed and 
there is suddenly a phase transition. 

Anticipating regime shifts in gene regulatory system 
(aka in protein concentration level) can be useful to pre¬ 
vent disease onset and progression which may intercept 
unacceptable sudden transitions from a healthy state to 
a disease state [l^. Examples of such regime shifts are 
asthma attacks, epileptic seizures and sudden deteriora¬ 
tion of complex diseases [l3, [13, [13| ■ A well documented 
example of regime shift is type 1 diabetes (TID) which 
is a form of diabetes mellitus [11]. TID is a chronic in¬ 
flammatory disease caused by insufficient production of 
insulin by j3 cells in the pancreas. The genetic association 
of TID is that the production of insulin through /3 cells 
depends on HLA-encoding genes M- If TID associated 
genes are in “on” expression state, then /3 cells produce 
insulin and releases insulin into the blood stream, but if 
they are in “off” expression state /3 cells fail to produce 
insulin which leads to TID. Prediction of regime shifts 


in gene expression using EWS from normal (“on”) to di¬ 
abetic (“off”) state could prevent the switch to diabetic 
state and help to maintain the level of insulin. 

In summary, our work reveals that stochasticity can 
have diverse complex effects on genetic regulatory sys¬ 
tems. Early warning signals to anticipate forthcoming 
regime shifts in gene expression requires special atten¬ 
tion to the underlying various statistical issues and limi¬ 
tations. In addition, to select a suitable window size and 
data length raise further difficulties. The main challenge 
of detecting early warning signals includes risk of false 
alarms and failed detections. One needs a deeper under¬ 
standing of early warning signals of regime shifts, and 
how a balance between early warning signals and false 
alarms is achieved, will lead to important new insights 
in genetic regulation. Moreover, our results establish the 
important fact that finding a more robust indicator of 
regime shifts in complex natural systems is still in its in¬ 
fancy and demands extensive research. We hope that this 
study may also increase the interest among researchers 
to find a more robust indicator for detecting upcoming 
sudden transition in much broader class of systems in 
developmental biology. 

SUPPORTING INFORMATION 

Appendix SI Derivation of Stationary Proba¬ 
bility Density Function. (PDF) 

Appendix S2 Additional examples of early warn¬ 
ing signals. (PDF) 
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